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Abstract 

Background: The Hill function and the related Hill model are used frequently to 
study processes in the living cell. There are very few studies investigating the 
situations in which the model can be safely used. For example, it has been shown, at 
the mean field level, that the dose response curve obtained from a Hill model agrees 
well with the dose response curves obtained from a more complicated Adair-Klotz 
model, provided that the parameters of the Adair-Klotz model describe strongly 
cooperative binding. However, it has not been established whether such findings 
can be extended to other properties and non-mean field (stochastic) versions of the 
same, or other, models. 

Results: In this work a rather generic quantitative framework for approaching such a 
problem is suggested. The main idea is to focus on comparing the particle number 
distribution functions for Hill's and Adair-Klotz's models instead of investigating a 
particular property (e.g. the dose response curve). The approach is valid for any 
model that can be mathematically related to the Hill model. The Adair-Klotz model is 
used to illustrate the technique. One main and two auxiliary similarity measures were 
introduced to compare the distributions in a quantitative way. Both time dependent 
and the equilibrium properties of the similarity measures were studied. 

Conclusions: A strongly cooperative Adair-Klotz model can be replaced by a suitable 
Hill model in such a way that any property computed from the two models, even 
the one describing stochastic features, is approximately the same. The quantitative 
analysis showed that boundaries of the regions in the parameter space where the 
models behave in the same way exhibit a rather rich structure. 



Background 

The Hill function and the related Hill model [1] are used frequently to study biochem- 
ical processes in the living cell In strict chemical terms Hill's model is defined as 

C + hA±*C h (1) 

where C denotes a protein that binds ligands, A is a ligand, and is a ligand-pro- 
tein complex having hA molecules attached to C. The stoichiometric coefficient h 
describes the number of ligand binding sites on the protein. All ligands bind at once. 
Both the forward and the back reactions are allowed. It is relatively simple to derive 
the expression for the dose response curve (the Hill function) which relates the 
amount of free ligands, a, to the fraction of ligand-bound proteins (e.g. receptors) in 
the system, 0. The Hill function is given by 
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a' 



cp = 




(2) 



where K 0 denotes the dissociation constant. 

The Hill function is used frequently in various areas of physics, biology, and chemis- 
try. For example, it is widely used in pharmacological modeling [2], as well as in the 
modeling of biochemical networks [3]. In the most common scenario, the Hill function 
is fitted to an experimentally obtained dose response curve to infer the value of the 
stoichiometry coefficient, h. The value obtained in such a way is not necessarily an 
integer number and is referred to as the Hill coefficient. The number of ligand binding 
sites is an upper limit for the Hill coefficient. The Hill coefficient would reach this 
limit only in the case of very strong cooperativity. More discussions on the topic can 
be found in [4]. However, in present study, the variable h will be allowed only non- 
negative integer values. 

Hill's model has been heavily criticized since it describes a situation where all ligands 
bind in one step [5]. In reality, simultaneous binding of many ligands is a very unlikely 
event. A series of alternative models have been suggested where such assumption is 
not implicit [6-8]. A typical example is the Adair-Klotz model [6] defined as 



with i = 1, h\ Protein C binds ligands successively in K steps. Here, and in the fol- 
lowing, the subscript i on C denotes the number of A molecules attached to it, with 
the obvious definition C 0 = C. Apparently, in comparison to the Hill model, the alter- 
native models - while being more realistic - are more complicated and harder to deal 
with (e.g. the Adair-Klotz model shown above). Accordingly, the central question being 
addressed in this work is whether it is possible to establish conditions where Hill's 
model can be used safely as a substitute for a more complicated reaction model. With 
a generic understanding of when this can be done, it should be possible to study an 
arbitrary reaction system with the elegance that comes with the use of Hill's model, 
knowing at the same time that the results are accurate. Also, even if there is evidence 
that the Hill model might describe the problem, it is not immediately clear which fea- 
tures of the problem can be described faithfully. 

In the following, Hill's model will be compared with a well chosen reaction model 
that is more realistic, and not too complicated from the technical point of view. The 
Adair-Klotz model discussed previously is a natural choice since it assumes that 
ligands bind sequentially, and the model is relatively simple to deal with. 

Furthermore, it is necessary to choose which property to study. For example, Hill's 
and Adair-Klotz's models have been compared in [5] where the property of interest 
was the dose-response curve (j){a). Using classical chemical kinetics, the dose-response 
curves predicted from Adair-Klotz's and Hill's model were compared neglecting fluc- 
tuations in particle numbers. It was found that for a strongly cooperative Adair-Klotz 
model it is possible to find the parameters for Hill's model that will result in similar 



Q_i +A 4Q 



(3) 



Pi 

Q -> Q_i+A 



(4) 
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dose response curves. The question is what happens for other properties, and what 
happens when fluctuations in particle numbers are taken into account? 

To avoid dealing with a particular choice of a property of interest, and to strive for 
an exact treatment, the models will be compared on the level of the respective particle 
number distributions. The position developed in this work is that the particle number 
distribution function of a model is the fundamental quantity that describes all features 
of the system. If the particle number distributions are similar, any property computed 
from them should have numerical values that are close. For example, the relevant vari- 
able for both models is the number of free ligands in the system. If the particle num- 
ber distribution functions are same for both models then the resulting number of free 
ligands will be same. However, the opposite might not hold: it might be that the num- 
ber of free ligands is same but some other quantity (e.g. fluctuations in the number of 
free ligands) might be be vastly different. To avoid such traps, the focus is on compar- 
ing the particle number distribution functions directly. 

The scope of the analysis in [5] will be extended in several ways. First, in addition to 
studying the stationary (equilibrium) properties of the models, dynamics will be studied 
as well. Many processes in the cell are strongly time dependent and involve coopera- 
tive binding, such as the early stages of signalling processes, and cascades in later 
stages of signal propagation phase. Likewise, many processes in the cell need to happen 
in a particular order. Clearly, the time and dynamics play a crucial role in the workings 
of cell biochemistry. Second, the previous mean field (classical kinetics) analysis will be 
extended to account for effects of fluctuations (intrinsic noise). It has been recognized 
that intrinsic noise (fluctuations in the numbers of particles) is not just a nuisance that 
the cell has to deal with, but is an important mechanism used by the cell to function 
[9-12]. Intrinsic noise becomes important when protein copy numbers are low. Such a 
situation is frequent in the cell (e.g. gene expression networks). Third, a generic com- 
parison of the models will be provided by focussing on the particle number distribu- 
tion functions. 

Results and discussion 

Description of models 

The models are parameterized as follows. Hill's model is parameterized by two reaction 
rates for the forward and the back reactions that will be denoted by a and /? respec- 
tively. The dissociation constant for the model K 0 is governed by the ratio p/a and for 
simplicity it will be assumed that 



The Adair-Klotz model involves more parameters: the forward and the back reaction 
rates for an i-th reaction are given by a t and p t respectively, and i = 1, h\ The disso- 
ciation constants for the Adair-Klotz model are defined as 



K 0 = P/< 



(5) 



Ki = -;i= 1, ft' 



(6) 



It is assumed that the particles mix well and that it is sufficient to count the parti- 
cles. The models are stochastic and are described using the continuous time Markov 
chain formalism [13]. The reaction rates govern the transition probabilities between 
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states of the system. The master equations for the models are the consequence of the 
corresponding forward Chapman-Kolmogorov equations for the transition probabilities. 
The solutions of the master equations are the particle number distribution functions as 
explained in the "Computation of the distribution functions" section. To compare the 
distribution functions for the models, three similarity measures are defined in the 
"Comparison of the distribution functions" and "Fine tuning the comparison proce- 
dure" sections. 

From the model-centric view taken in this investigation, the best way to compare the 
distribution functions is to choose h - h\ This makes the number of binding steps in 
the Adair-Klotz model equal to the stoichiometric coefficient of the Hill model. Also, 
within the scope of this work, to simplify wording, the variable h will be simply 
referred to as the Hill coefficient. The choice h = K makes it possible to relate the dis- 
tribution functions in a rather natural way. Namely, if h = h\ it is possible to establish 
a one to one correspondence between Hill's model state space and a subspace of 
Adair-Klotz's model state space. The respective states in these spaces will be referred 
to as common states, or the common state space. 

The first similarity measure defined, S(t), quantifies the similarity between the distri- 
bution functions for Hill's and Adair-Klotz's models on the space of common states. In 
the text this similarity measure is referred to as the main or fundamental similarity 
measure. The states in Adair-Klotz's model state space that are not part of the com- 
mon state space are referred to as the complement (state) space. This set contains 
states in which at least one of the intermediate species (section "Computation of the 
distribution functions") is present. These states are unique to Adair-Klotz's model. 

The second similarity measure introduced <5(t), measures the extent to which the 
complement space is occupied. This is an auxiliary similarity measure that comple- 
ments the information conveyed by the use of the fundamental similarity measure d{t). 

The third similarity measure, $ (t),) quantifies the similarity between the shapes of 
Hill's and Adair-Klotz's model distribution functions. It is also an auxiliary similarity 
measure used to refine the information provided by inspection of the fundamental 
similarity measure. To compare the shapes of the distribution functions, Adair-Klotz's 
model distribution function is re-normalized on the common state space. 

Optimization of Hill's model parameters 

One needs to be careful not to compare an arbitrary Hill's model to an arbitrary Adair- 
Klotz's model. Since the goal is to quantify which Adair-Klotz's models can be replaced 
by the related Hill's models, it is natural to choose the best possible parameters for the 
Hill model that maximize the fundamental similarity measure S{t). Thus for each 
choice of the parameters for the Adair-Klotz model, the parameters of the Hill model 
will be optimized. The optimization procedure differs somewhat for plots that depict 
time dependence from the ones that depict equilibrium properties. 

In the equilibrium, S(t) depends only on the values of the dissociation constants: 
- lim^ >OQ S(t) and 

8 oo =fiK 0 ,K 1 ,K 2 ,...,K h ) (7) 

For a fixed tuple (K lf K 2 > K^) the Hill model dissociation constant K 0 is optimized 
to make cL as large as possible. This makes the Hill's model dissociation constant 
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dependent on Adair-Klotz's model dissociation constants in a well defined way: 

K 0 =g{K lf K 2 ,...,K h ) (8) 

where g is the function resulting from the optimization procedure. Thus one can 
write ^ = f(g(Ki, K 2 > Kh)> Kit — > Kh)> which defines the function ^ max such that 
for a given choice of dissociation constants for the Adair-Klotz model 5^ is the largest 
possible 

Soo = &mB X (K 1 ,K 2 ,...,K h ) (9) 

The function <5 max is depicted in all plots that analyze the equilibrium state. 

Please note that the use of Eq. (8) only fixes the ratio p/a. Accordingly, for time 
dependent plots, an additional choice has to be made for either a or /5. For a time 
dependent plot the value for a was adjusted so as to make the life-time of the initial 
state the same in both models. (During the optimization, the value of /3 is given by 
K 0 a). 

Numerical results 

The three similarity measures have been computed numerically by solving the master 
equations for the models. Figure 1 shows how the similarity measures 
A(t) e {8{t) f 8(t), 8{t)} depend on time in the situation where it is expected that Hill's 
model cannot approximate the dynamics of Adair-Klotz's model, i.e. when all reaction 
rates are equal and Adair-Klotz's reaction system cannot be described as cooperative. 
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Figure 1 Similarity measures (weakly cooperative Adair-Klotz model, h = 2). Time dependence of the 
similarity measures for h = 2 case: A(t) G {5(f), 8{t), 8{t)} This and all other figures in the 
manuscript were generated with P 0 = 2 and L 0 = 5. In this figure weakly cooperative Adair-Klotz model has 
been considered with a, = ft, = 1s" 1 for / = 1, h. The parameters for the Hill model were optimized so 
that S(t) is largest possible {p/a = 0.5 and a = 0.55" 1 ). The time f is expressed in units of s. The full line is 
for A = S, while the dashed and the dotted lines are for ^ = § and A = 5 respectively. 
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The similarity is perfect at t = 0 by construction, since in principle both systems are 
prepared in identical states. The similarity starts decreasing since the intermediate 
states become populated. This can be seen from the fact that the dashed line goes up, 
starting from zero. Please note that after some time the intermediate states become 
de-populated since the dashed line goes down after the initial peak around t « 0.25. 
The choice of reaction rates for the Adair-Klotz model clearly makes the intermediate 
states long lived. In such a case it is not possible to find the parameters a and /? such 
that the fundamental (main) similarity measure is large. 

The first auxiliary similarity measure that relates the shapes of the distribution func- 
tions (the dotted line in the figure) exhibits interesting behaviour: §(t) & 1 for all times 
(early, intermediate, and asymptotic). Given this insight, one can conclude that only 
properties (observables) that are shape sensitive can be described by Hill's model, 
despite the fact that intermediate states are highly populated. For example, the 
moments of the particle number distributions do not fall into this category (e.g. the 
average numbers of particles in the systems or the variances); however, ratios of 
moments (defined on the common state space) do. 

To which extent are the findings discussed so far sensitive to the value of the Hill 
coefficient? Figure 2 was constructed in the same way as Figure 1, but with a higher 
value of the Hill coefficient. To make the computations faster, the lowest possible 
value for the Hill coefficient was used, i.e. h = 3. In comparison to the h = 2 case, the 
fundamental similarity measure decreases further. It can be seen that 8(t) increases, 
which indicates that the complement space becomes more populated. It is very likely 
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Figure 2 Similarity measures (weakly cooperative Adair-Klotz model, h = 3). Generated in the same 
way as Figure 1, but with a higher value for the Hill coefficient (h = 3). The parameters of the Hill model 
were optimized in the same way as for Figure 1, resulting in a = 0.5s" 1 , /3 = 0.083s" 1 . Increase in the Hill 
coefficient makes the discrepancy larger since there are more intermediate states that can be populated. 
The similarity in the distributions shape increases for large times. 



Konkoli Theoretical Biology and Medical Modelling 201 1, 8:10 
http://www.tbiomed.eom/content/8/1 /1 0 



Page 7 of 1 7 



that this is because more intermediate states are available. The shape similarity mea- 
sure 8 (t) decreases for intermediate times, as the dotted curve has a deeper minimum 
than the dotted curve in Figure 1. 

For the case in which intermediate states are short lived, one intuitively expects that 
Hill's model could be a useful substitute for Adair-Klotz's model. Figure 3 depicts the 
dependence of the similarity measures on time, for systems that are expected to behave 
in a similar way. In particular, the reaction rates for the Adair-Klotz model used were 
chosen in such a way that the intermediate states are short lived. Indeed, the value of 
8{t) stays very close to 0. The shapes similarity measure 8 (f) stays very close to one, 
finally leading to large values for the fundamental similarity measure 3{t). This is an 
important finding since it indicates that Hill's model can be used to investigate an arbi- 
trary observable, e.g., not just the average number of free ligands, but also the noise 
characteristics of that quantity. Naturally, such a claim comes with the implicit con- 
straint that the observable should be interpreted in the context of Hill's model state 
space. For example, quantities such as the number of free receptor proteins, or the 
number of fully occupied receptors, fall in this category. However, any quantity that 
would involve counting the number of intermediates does not. 

The time dependence of the similarity measures was investigated to confirm that 
these analysis tools work as expected. It is important to check that the analysis will 
work for both dynamics and the equilibrium state. In the following, the focus is on 
understanding equilibrium properties. The goal is systematically to identify situations 
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Figure 3 Similarity measures (strongly cooperative Adair-Klotz model, h - 2). Generated in the same 
way as Figure 1, but with different values for the reaction rates. The particular choice of the reaction rates 
makes the intermediate states weakly populated: a } = 1s" 1 , fa = 10s" 1 , a 2 = 10s" 1 , and fa = 1s" 1 . The 
parameters for the Hill model were optimized in the same way as for the Figure 1 resulting in a = 0.5s" 1 
and ft = 0.25s" 1 . S{t) stays relatively close to one indicating a good match. The dashed curve stays low, 
which indicates that intermediate states are short lived. The dotted line stays close to one indicating that 
the distributions have a similar shape. 
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when Hill's and Adair-Klotz's model distribution functions are similar. Technically, this 
will be done by mapping out regions in the Adair-Klot's model parameter space where 
the fundamental similarity measure ^ max is relatively high. 

Figure 4 shows how ^ max depends on the values of the Adair-Klotz model reaction 
rates for the case h = 2. The figure depicts contours where S max = const in the (K lt K 2 ) 
plane. The first interesting region is in the range 0 < K\ ^ 45 and below the full curve. 
In this range (the grey region below the full curve) K\ » K 2 guarantees high similarity 
measure values. This analysis confirms the previous mean field study [5] where it was 
shown that choosing K x » K 2 leads to similar dose response curves. In the present 
article it has been shown that the results holds for any observable (average numbers, 
variances, etc). The second interesting region is for K x £ 45. In that region the funda- 
mental similarity measure is large for any K 2 . Cases with relatively large values of K 2 
are not interesting chemically, since such reactions would be chemically non-func- 
tional: K±K 2 » 1 would lead to the situation where the fraction of final products 




Figure 4 Equilibrium state similarity measure for h = 2 The contour plot that depicts how long time 
limit of <L = lim f _>oo S{t) depends on the dissociation constants K } = /V#i and K 2 = P2/&2, = f{K 0l K h 
K 2 ). For a fixed pair {K h K 2 ) the Hill model dissociation constant K 0 = fi/a is optimized to make <L as large 
as possible, making the Hill's model dissociation constant dependent on Adair-Klotz's model dissociation 
constants in a well defined way; K 0 = g(K h K 2 ) leading to the function = f(g(K h K 2 ), K h K 2 ) = S max {K h K 2 ) 
that is depicted in the plot. 
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(complexes) in the system would be vanishingly small. However, a reaction with K± £ 
45 and K 2 K< 1 could be functional provided KiK 2 ~1. 

Figure 5 shows similar kind of analysis as done for Figure 4 but for the first higher 
value of the Hill coefficient, h = 3. Unfortunately, because the structure of the para- 
meter space is more complicated, it is not possible to use a single contour plot. 
Instead, various hyperplanes in the parameter space are studied. Panel (a) depicts the 
regions in the (K lf K 2 ) plane where ^ max = 0.9 for different choices of 7<T 3 . The region 
with S max > 0.9 is always to the right of each curve. For example, in the grey region in 
panel (a), for K 3 = 1000, it is always true that ^ max > 0.9. On the one hand, it can be 
seen that increase in K 3 reduces the area where the fundamental similarity measure is 
large. On the other hand, for a fixed value of K 3 , and for a chemically functioning reac- 
tions {I<iK 2 choosing K x » K 2 makes the fundamental similarity measure large. 
Likewise, panel (b) indicates that to obtain a large value for the fundamental similarity 
measure Ki should be as large as possible. For a given value of Ki one should take K 2 
» K 3 . In brief, one can say that K\ » K 2 » K 3 ensures that ^ max is large but the plot 
shows that there are many subtle details associated with such a statement. Again, this 
confirms the previous finding in [5] that Ki » K 2 » K 3 results in similar dose 
response curves for both models, but please note that the statement made in here is 
much more general. 

The quantitative analysis reveals rather rich structure of the parameter space where 
the two models have very similar noise characteristics (distribution functions). It would 
be useful to simplify such criteria. In that respect, it is tempting to express the strong- 
cooperativity criteria 



in another way, e.g. by introducing a measure of the degree of cooperativity d? as 



The strong cooperativity can be characterized by df » 1. Naively, one would expect 
that in such a way one should obtain high values for S mSLX uniformly in K lt 

Figure 6 is a contour plot that depicts how ^ max depends on K ± and df for h = 4. The 
figure shows that many parameter choices that are chemically interesting do lead to a 
high value of the fundamental similarity measure (the grey region in the plot). Since 
there is no upper limit for df, for any value of K lf it is possible to choose df so that the 
reaction is chemically operational: for large d? the product KiK 2 K^K^ ~ K\/% 6 becomes 
very small. However, there is rather large region close to the origin (the white region 
in the plot) where the Hill model is not a good replacement for the Adair-Klotz 
model. The minimal value of df that guarantees a good match needs to be adjusted 
depending on a value of K lt Interestingly, for K x £ 65 any value of df will lead to large 
^max- Unfortunately, it was not possible to generate similar figures for h > 5 owing to 
the limitations of the computer hardware. 

Conclusions 

Particle number fluctuations as predicted by Hill's and Adair-Klotz's model have been 
studied quantitatively. To compare the fluctuation characteristics of the two models, 



Ki » K 2 » . . . » K h 



(10) 



(Ki,K2 K h ) = (K lf -± 



(11) 
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Figure 5 Equilibrium state similarity measure for h = 3 The plot depicts equilibrium state similarity 
measure for h = 3 case. For each triple {K h K 2l K 3 ) an optimal value is found for K 0 that maximizes <L. In 
such a way <L = d msK (K h K 2l K 3 ). The lines plotted in both panels denote the cL = 0.9 boundaries. For a 
given curve, the region with <L> > 0.9 is always to the right of the curve. Panel (a): the reaction rates 
parameter space is projected on to {K h K 2 ) plane with K 3 fixed at the values indicated in the panel. Panel 
(b): the parameter space is projected on the {K 2l K 3 ) plane with several choices for K } as indicated in the 
panel. 
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Figure 6 Validity region of a » K 2 » K 3 » K 4 parameterization The plots depicts the boundary of 
the S max {K h K 2l K 3l K 4 ) > 0.9 region in {K h Q plane with the parameterization K 2 = K-\/£ K 3 = K^^ 2 , and K 4 = 
^/(f 3 . {K 0 has been optimized as in the previous figures.) 



the similarity between the particle number distribution functions was characterized by 
three quantitative measures of similarity. The fundamental similarity measure S(t) 
expresses the degree of overlap between the distribution functions on the common 
state space. Two auxiliary similarity measures 8[t) and 8 (t) have been introduced to 
refine the analysis further by measuring the degree of occupancy of intermediate states, 
and measuring the similarity in the shape of the distributions on the common set of 
states. 

It was shown that the similarity measures work as expected by studying their time 
dependence. The value of S(t) always follows 1 — 8 (t)- This quantifies the intuitive 
expectation that the occupancy of the intermediate states governs whether models 
behave in the same way. In addition, it was found that, interestingly, 8 (t) stayed rela- 
tively close to one, even when S(t) was relatively small. 

Furthermore, the equilibrium similarity measure = lim^oo S(t) was analyzed, 
where dependence of 5^ on values of the dissociation constants K lf K 2 , K h was care- 
fully investigated. The analysis revealed that a value of the similarity measure in the 
equilibrium state is high when K x » K 2 » ... >:> K h . This is in agreement with findings 
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in an earlier work [5], which showed that the dose response curves for both models 
agree in this regime, provided the condition on the dissociation constants holds. 

This work extends previous findings by avoiding the mean field approximation, 
and focussing on the distribution functions. By doing so it is possible to extend the 
previous finding to any property of interest that can be obtained from the particle 
number distribution functions. Furthermore, it was shown that the boundaries of the 
parameter space where is high have a rather rich structure. While it is true that 
the condition I<i» K 2 » ... K h guarantees that a given Adair-Klotz model can be 
substituted by a Hill's model, there are subtle details that need to be attached to 
such a statement. 

The findings of this work should shed some light on the applicability of the previous 
uses of Hill's model. For example, Hill-like models have been used in the past to study 
characteristics of fluctuations in particle numbers during the process of complex for- 
mation [14-16]. This study shows that findings in these studies can be extrapolated to 
more realistic reaction models of complex formation, without doing the advanced tech- 
nical analysis required for understanding more realistic reaction models. 

This work can be extended in many ways. First, it should be possible to consider 
more challenging limits, with larger values of the Hill coefficient and particle copy 
numbers. Relatively small values for these parameters were considered owing to the 
limitations of the computer hardware (memory and CPU). Likewise, only pure states 
were considered, and it would be interesting to see whether the same conclusions can 
be drawn for other types of initial conditions. Second, instead of analyzing the full dis- 
tribution functions, it should be possible to investigate the similarity of the underlying 
moments, and to define similarity measures accordingly. This could be advantageous 
for studying the problematic limits discussed above. Third, the similarity with, and 
among, other reaction models could be studied in a way similar to that presented here. 
For example, the issue of model reduction is a perpetual everlasting problem in the 
modelling of intracellular processes. 

Methods 

Computation of the distribution functions 

To compare the models the particle number distribution functions will be investigated. 
It will be assumed that particles mix well. In such a setup, it is sufficient to count the 
particles. The numbers of C 0 , C lf C 2 , C h and A particles will be denoted by n 0 , 
n 2 , rih and n A respectively. 

Each system has a configuration space associated with it. The configuration spaces of 
the system are similar but not identical. For Hill's model a configuration of the system 
is given by c H = (n 0 , n h , n A ), while for Adair-Klotz's model c A = (n 0 , n lf n 2 , n h) n A ). 
The difference comes from the fact that molecules Q, C 2 , C h _ x need to be counted. 
In the following these molecules will be referred to as the intermediate molecules or, 
in brief, the intermediates. 

The systems are stochastic and in course of time transitions within the configuration 
spaces of the systems occur randomly. The rapidity of transitions is governed by the 
previously introduced reaction rates. Both systems can be described by their respective 
master equations. 
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The master equation for Hill's model is given by 

Wch, t) = a K + h ) P H (c H [+, -, +], t) 

+P(n h+ l)P H (c H [-, +l -],t) (12) 



■(?) 



where d t denotes the time derivative. The states c H [+,-, +] and c H [-,+r] are defined 
by. 

c H [±, ±, ±] = (n 0 ±l,n h ± 1, n A ± ft) (13) 

where any combination of the plus and the minus signs can be picked at will (a 
choice has be to made consistently by picking either all upper or all lower signs). The 
particle number distribution function Ph{ch> t) defines the occupancy probability for a 
state c H at a time t. 

The master equation for Adair-Klotz's model is given by 

h 

d t P A (c A , t) = ^2 K(n f _i + l)(n A + l)x 

1=1 (14) 

xP A {c A [i, +], t) + Pi{m + i)P A (c A [i, -], t) 
-{ptini-iriA + Pirii)P A {c A , t)] 

where 

c A [i, ±] = (n 0 , • • • , ± 1, n f =f 1, . . . , n A ± 1) (15) 

where either the upper or the lower set of signs can be picked at will. 

By solving the master equations (12) and (14) it is possible to obtain the distribution 
functions P H and P A for Hill's and Adair-Klotz's models respectively. In the next sub- 
section the procedure for comparing the distributions will be discussed. 

Structure of the configuration spaces 

To make a fair comparison between the models it is natural to use the same initial 
conditions for both. Since Hill's model does not have information about the intermedi- 
ates, the initial conditions will be chosen so that the copy numbers of the intermediate 
species are all zero. 

For Hill's model the dynamics will be started from a pure state with initial configura- 
tion given by 

c° H = {Po,0,L 0 ) (16) 

where P 0 and L 0 denote the number of protein complexes and the number of ligand 
molecules in the system at t = 0. Likewise, for Adair-Klotz's model, the system will be 
started from 

c° = (Po,0,...,0,L 0 ) (17) 

For the pure initial state the dynamics of Hill's model occurs on the one dimensional 
space defined by the following states 
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4 = ( p o-i,U 0 -ta) ( 18 ) 

where i = 0,1,2, . . . , i 1 ^ and the upper limit for the state index i is given by 
i H = mm(L 0 /h, Po). The initial state corresponds to i = 0. This set of states will be 
referred to as 

S„ = {4|i = 0,l,2 Cox) ( 19 ) 

Likewise, for a pure initial state, following set of states emerge for the Adair-Klotz 
model, 

cY 2 lh = [Po ~ [h + ii + • • • + ih), 

ii#i2/ • • • / ih/ ( 2 °) 
L 0 - (z'i + 2i 2 + . . . + Wh)) 

Such set of states will be referred to as the Adair-Klotz space and denoted by 

Sa = {cY 2 fh |ii,f2,-.-,ih = 0 / l,...*} (21) 

where symbol * in the equation indicates that the upper limit has to be chosen such 
that occupancy numbers for each configuration are positive. Equation (20) indicates 
that protein molecules are either free from ligands, or have one or more ligands 
attached to them. From the perspective of the ligands, the equation states that all 
ligands that are not free are bound to protein molecules either as a single molecule, or 
in pairs, triples etc. 

The inspection of the configurations for Hill's and Adair-KlotzVs models, in (18) and 
(20), reveals that the configuration spaces are rather similar, up to the fact that the 
Adair-Klotz space has much higher rank. 

Furthermore, it is possible to see that a vector in Adair-Klotz space with i x - 0, i 2 - 
0, 4_! = 0 (Eq. 20) has a natural correspondence with the vector in the Hill space 
with i = i h (Eq. 18). In what follows it will be useful to formalize this mapping. 

Symbol £a(ch) will denote the image of a state c H in the Adair-Klotz space, 

l A (c H ) = (n 0 , 0, 0, . . . , 0, n h , n A ) (22) 

The set of images of all vectors in the Hill space will be denoted by 

X A {S H ) = {X{c H )\c H e S H } (23) 

Please note that this mapping defines a one to one correspondence between the 
states in the Hill and the Adair-Klotz spaces. For example, given that i and h are fixed, 
there is only one combination of i lt ... i h for which i = i x + i 2 + ...+ ih and hi = i\ + 2i 2 
...+ hi h . 

Clearly, X A (Sh) C S a , and the set of states that are in the Adair-Klotz space but not 
in the image space (i.e. a complement) will be denoted by C A (Sh) = S A \T A (Sh). 

Comparison of the distribution functions 

To compare the probability distributions for the models, the distribution function for 
Adair-Klotz's model will be projected on to the state space of Hill's model: 

Pa{c h , t) = Pa{X a {c h ), t) (24) 
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The direct comparison of P H (c H ) with P a (ch) can reveal whether there is a region in 
the parameter spaces of the two models where the respective dynamical behaviour is 
similar. 

Once the projection is done, the comparison of the distribution functions is equiva- 
lent to the comparison of two vectors in a Cartesian space. For example, it is possible 
to use the scalar product between the vectors to compare them. However, for the pur- 
pose of this work, the distributions will be compared using 



8{t) = ^ p h{c Hi t)P A {c Hl t) 



(25) 

c H eS H 

The advantage of the particular form used in (25) is that for the perfect match with 
Ph{ch) = Pa{ch) for all c H e S H , the similarity measure S(t) equals one. This can be 
seen from that fact that the sum in (25) becomes the normalization condition for the 
distribution functions. The lowest value for S(t) is clearly zero since the distribution 
functions are positive definite. Also, please note that in the light of (16) and (17), ^(0) 
= 1. The initial conditions are chosen so that the match is perfect at t = 0. In such a 
way, any discrepancy detected by S(t) is due to the dynamics of the systems. 

Fine tuning the comparison procedure 

In addition to the similarity measure defined in Eq. (25) it is useful to analyze the 
extent to which the states in the complementary space Ca (Sh) are populated. In that 
respect, it is useful to introduce 

*W = P ^> 0 (26) 

c A eC A {S H ) 

This measure is important since it indicates to what extent the presence of inter- 
mediates affects the value of S(t) in (25). 

If the intermediate states are short lived, they should not be populated, and accord- 
ingly 8(t) ~ 0. In such a case S(t) has a fair chance of being equal to one. On the other 
hand, for 8(t) ~ L S(t) will be small, although the fact that the shapes of Hill's model 
distribution and Adair-Klotz's model distribution (projected on S H space) might be 
similar. 

To analyze quantitatively the effects discussed above, it is useful to introduce a mea- 
sure of the similarity of Hill's model distribution function and the normalized distribu- 
tion function of Adair-Klotz's model P{c H , t) on Hill's space. To do this, it is useful to 
renormalize Adair-Klotz's model distribution function on the image space as 

D ( +\ ^a(ch/0 

p a{ch, t) = jz jr (27) 

\PA{c H ,t)\ 

where the norm is given by 

\h(c H , 0 1| = J2 Pa ( Ch ' f ) (28) 

c H eS H 



Please note that since Adair-Klotz's model distribution function is normalized, the 
following condition holds 
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\Pa{c h , t)|+5(t) = l (29) 

The similarity measure of Hill's model distribution function P H and the renormalized 
distribution function of Adair-Klotz's model p A can be finally defined as 

Ht) = J2 \/ p h(ch, t)P A {c H , t) (30) 

Please note that 8 (t) measures the similarity in the shapes of the distribution func- 
tions constrained on the Hill space, and in this work is referred to as the shape simi- 
larity measure. 

Finally, using the equations above, it is trivial to show that 

8{t) = yjl - 8(t)8(t) (31) 

The similarity of distributions can be factored in two contributions. The square root 
term on the right hand side of the equation measures the extent to which the image of 
the Hill space is populated for Adair-Klotz's model. The second term on the right 
hand side of the equation measures the similarity of the shape of the probability distri- 
butions on Hill's space image. To obtain a good match, both factors in the product 
need to be large, the intermediates should be short lived, and the shape of the distribu- 
tions should be similar. 



Numerical computation setup 

The distribution functions were computed by Mathematica using the technique of the 
Laplace transform. The Laplace transform of a function fit) is defined in the usual way 
as 

poo 

C [f{t) t s] = / dte-* f[t) (32) 
Jo 

The Laplace transform of the time derivative becomes an algebraic expression. Using 
this property, a master equation can be converted into an algebraic equation. The 
resulting linear algebraic equations were solved using the internal Mathematica solver. 
The asymptotic time limits of time-dependent functions were computed easily using 

lim/(t) = lim sC[f{t), s] (33) 

t^OO 5^0 

Accordingly, the equilibrium quantities were computed with infinite precision. 

For the time dependent quantities, the numerical inversion of the Laplace transform 
for the distribution functions was done using the Durbin method. The computations 
were performed using the Mathematica package developed by Arnaud Mallet and can 
be found at the repository of Mathematica packages. Thus the numerical results shown 
in the figures for time dependent quantities are exact to the accuracy of the numerical 
inversion procedure. The inversion formula is based on an integral that needs to be 
evaluated numerically. The accuracy of the result depends on the number of points 
used to perform the integral. This number was doubled incrementally until the relative 
change in the computed value was below 1%. 
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